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I.  INTRODUCTION 


Laser  guidance  systems  operating  in  the  atmosphere  are  affected 
by  several  phenomena  which  tend  to  degrade  the  transmitted  laser  beam 
and  result  in  poorer  performance  than  would  be  obtainable  in  vacuum 
transmission.  Of  these  phenomena,  the  major  contributors  to  beam  degra- 
dation are  direct  attenuation  of  the  beam's  energy  by  absorption,  spread-  i 

ing  and  motion  of  the  beam  due  Co  random  index-of-refraction  variations 
within  the  atmosphere,  and  scattering  of  the  electromagnetic  radiation 
by  atmospheric  particles.  i 

Direct  attenuation  of  the  transmitted  energy  occurs  from  molecular  I 

absorption  by  the  various  constituents  in  the  atmosphere  and  absorption 
by  aerosol  particles  present  in  the  atmosphere.  Molecular  absorption 
occurs  at  selected  lines  or  frequencies  and  varies  greatly  over  fairly 
narrow  frequency  intervals.  Aerosol  absorption,  however,  remains  nearly 
constant  for  a given  particle  over  wide  ranges  in  frequency  of  the 
transmitted  energy.  Methods  are  presented  in  this  report  for  determining 
the  absorption  coefficients  for  molecules  and  aerosols. 

Index-of-refraction  variations  in  the  atmosphere  occur  from  heat 
transfer  processes  which  produce  air  temperature  Inhomogeneity  or  !' 

turbulence.  Movement  of  these  random  fluctuations  in  air  temperature 
across  the  path  of  a transmitted  laser  beam  caused  by  wind  or  beam  ! 

sluing  result  in  the  beam  wandering  and  spreading  about  its  original 
aim  point.  Thus,  the  energy  distribution  or  "spot"  produced  at  a dis- 
tance downrange  of  the  transmitter  does  not  remain  constant,  but  enlarges 
and  wanders  about  the  "target"  area.  This  turbulence- induced  beam  |; 

"jitter"  or  wandering  adds  an  additional  component  to  that  already  pro- 
duced  in  the  guidance  system.  This  report  describes  a new  analytical  j; 

method  for  rapidly  predicting  spot  movement  and  size  as  a function  of  |i 

transmitter  characteristics,  refractive  index  structure  constant,  and  i; 

effective  wind  speed  across  the  beam.  | 

Air  molecules  and  other  particles  such  as  dust,  haze,  fog,  or  smoke, 
which  are  present  in  the  atmosphere,  degrade  a laser  beam  by  scattering 
part  of  the  energy  out  of  the  path  of  the  beam.  For  particles  very 
small  relative  to  the  laser  wavelength,  the  angular  distribution  of  the 
scattered  energy  can  be  easily  determined.  For  particle  sizes  of  the 
same  order  as  the  wavelengtli  or  larger,  however,  the  theory  becomes  more 
complex  and  extensive  calculations  are  required.  Analytical  methods 
suitable  for  making  these  calculations  are  well  developed  for  most 
particles  of  Interest,  but  the  physical  data  required  are  not  adequate 
in  many  cases.  For  example,  the  complex  indlces-of-refraction  and 
particle  size  distributions  of  many  aerosols  are  not  well  known. 

This  effort  was  undertaken  to  provide  methods  that  couid  be 
used  to  make  predictions  of  the  effects  of  the  atmosphere  on  a terminal 
homing  laser  guidance  system.  Available  procedures  and  methods  were 
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reviewed  and  those  adequate  for  the  present  application  have  been 
adapted  for  use.  Scientific  support  was  obtained  to  provide  help  in 
developing  new  procedures  and  modifying  existing  models.  In  particular, 
the  efforts  of  Dr.  D.  L.  Fried  of  Optical  Sciences  Consultants  provided 
the  basis  for  the  turbulence  beam  wander  model  and  the  support  of 
W.  G.  Blattner,  D.  G.  Collins,  and  M.  B.  Wells  of  Radiation  Research 
Associates  was  obtained  to  modify  their  existing  Monte  Carlo  radiation 
transport  model.  The  molecular  and  aerosol  attenuation  models  were 
obtained  from  Dr.  A.  Miller  of  New  Mexico  State  University  and  Dr.  R.  B. 
Gomez  at  the  Army  Atmospheric  Sciences  Laboratory.  An  additional  model 
for  predicting  first-order  radiation  transport  for  a laser  designator 
system  was  developed  in-house. 

Descriptions  are  presented  in  the  following  sections  of  the  models 
developed  for  predicting  atmospheric  effects  encountered  by  a trans- 
mitted laser  beam,  and  the  Appendix  provides  utilization  instructions, 
Fortran  listings  of  the  in-house  developed  procedures  are  included. 


II.  TURBULENCE-INDUCED  BEAM  WANDER  MODEL 

For  predicting  turbulence- induced  wander  of  the  laser  trans- 
mitter beam  the  model  developed  by  Fried  [1]  for  calculating  the  power 
spectrum  of  angle- of-arrival  fluctuations  is  used  in  a numerical  pro- 
cedure employing  the  fast-Fourier  transform  (FFT)  to  convert  the  fre- 
quency dependent  power  spectrum  into  the  time  domain  giving  the  angular 
jitter  of  the  beam.  In  this  method,  the  power  spectrum  defined  as  the 
Fourier  transform  of  the  angle-of-arrival  temporal  coherence  function 
is  expressed  as 


F„(f)  - 4 y 

o 


cos(2jrfT)C^(T)dT 
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where  C^(t)  “ <Q:(t)Q!(t  + •r)>  is  the  temporal  coherence  function  of  the 
angle-of-arrival  fluctuations. 


In  the  numerical  procedure,  the  turbulent  region  between  the  trans- 
mitter and  target  (Figure  1)  is  subdivided  into  N segments  of  length 
and  the  calculation  performed  for  each  segment  is 


iJ 
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where 


\ * the  laser  wavelength 

» the  spot  diameter 

r^  * the  coherence  length  of  segment  i defined  by 


f^  to  a reference  frequency  given  by 


eff.l 

^o.i 71^ 

\ s t 


G^(f/fo  is  a function  which  can  be  approximated  by 


, if  0 S f < 0.332  f 


, if  0.332  f , < f < 3.1  f , 
o,i  o,i 

, if  3.1  f . < f 
o,i 


^ is  the  refractive  index  structure  constant  of  segment  i and  is 

the  distance  from  the  virtual  point  source  to  the  center  of  segment  i. 
^ is  the  effective  wind  velocity  across  the  beam  at  segment  i and 

can  include  the  effects  of  beam  sluing  by  combining  with  the  actual 
crosswind  velocity  according  to 

".tf,!  ■ ’..I  * - h'>  • 

where  the  plus  or  minus  sign  is  chosen  to  account  for  the  wind  being 
opposed  to  or  in  the  same  direction  as  the  angular  sluing.  The  power 
spectrum  for  frequency  f is 


- I 
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The  variance  of  the  power  spectrum  is  the  integral  over  frequency,  or 


M 

■ I ■ 

i-0 

where  M is  the  number  of  frequency  values  chosen. 

Values  obtained  for  the  power  spectnim  are  next  combined  with  a 
set  of  random  values  chosen  so  as  to  have  zero  mean  value  and  unity 
variance.  This  gives  a random  sequence  of  values  having  the  same  vari- 
ance as  the  calculated  power  spectrum.  The  relation  used  to  form  the 
random  sequence  is 

\<'l)  ■ “r,!  • 

where  N . is  one  of  the  normally  distributed  random  values  and  At  is 
^ * i 

the  time  interval  between  values  in  the  time  domain  chosen  so  that 
MM  At  - 1 

A symmetric  array  of  2M  values  is  obtained  by  folding  the  N^(fj^)  array 

to  give  an  equal  number  of  negative  frequency  values.  Fourier  trans- 
forming the  resulting  2M  array  gives  an  ordered  sequence  of  values 
representing  the  timj  Interval  -T  - 0 £ T,  where  T “ MAt.  The  second 
half  of  this  sequence  represents  one  component  of  beam  jitter.  A second 
set  of  independent  values  is  obtained  in  the  same  manner  for  the  other 
component  of  angular  jitter  to  give  spot  centroid  motion  at  the  target. 

To  obtain  beam  jitter  for  the  target- to- seeker  path,  the  power 
spectrum  of  angle- of-arrival  fluctuations  is  determined  by  a similar 
procedure  except  that  angle-of-arrlval  isoplanatism  effects  must  be 
considered.  The  complete  procedure  has  been  discussed  in  detail  by 
Fried  [2].  Total  power  spectrum  for  the  two  paths  becomes 

FT(f)  - F^(f)  + F^(f)  , 

where  F^(f)  is  the  power  spectrum  associated  with  the  viewing  process 
from  the  seeker. 

The  method  for  target  spot  size  determination  proposed  by  Fried  [2] 
is  to  make  use  of  his  short- exposure  resolution  theory  for  the  turbulence- 
induced  and  diffraction-limited  beamspread.  The  calculation  of  the 
effective  beamspread  becomes 
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where  6 “ 1.128  \/r  and  the  function  iI>(Dj/r  ) Is  determined  by  eval- 

r o u o 

uating  the  integral 


where  r^  is  the  coherence  length  for  the  designator  path. 

Total  angular  beamspread  is  determined  by  combining  with  the  trans- 
mitter optical  divergence  according  to 


and  the  corrected  spot  size  at  the  target  including  turbulence,  diffrac- 
tion, and  designator  optics  becomes 


Comparisons  of  the  results  of  calculations  of  beamspread  made  by 
this  method  with  the  methods  of  others  have  shown  differences  in  spot 
size  which  become  quite  large  under  conditions  of  strong  turbulence. 
Resolution  of  these  differences  to  determine  a preferred  method  of 
calculating  spot  size  awaits  an  experimental  validation. 


III.  MOLECULAR  AND  AEROSOL  ATTENUATION  MODELS 

Methods  for  the  prediction  of  molecular  and  aerosol  absorption 
and  scattering  developed  for  other  programs  were  adapted  for  use  in  the 
present  effort.  Models  developed  for  the  Army  Atmospheric  Sciences 
Laboratory  by  Miller  et  al.  [3]  and  Gomez  et  al.  [4]  are  particularly 
applicable  to  the  kinds  of  problems  considered.  Brief  descriptions  of 
these  two  models  are  presented  in  this  section.  Input  data  instructions 
for  making  calculations  are  Included  in  the  Appendix. 


li 


i;. 


i. 

i, 


8 


The  molecular  absorption  and  scattering  model  [3]  uses  the  AFGL 
(formerly  AFCRL)  line  parameters  compilation  [5j  and  computes  high 
resolution  molecular  absorption,  molecular  continuum  absorption,  and 
Rayleigh  scattering  coefficients.  Property  data  for  the  1962  US  Standard 
Atmosphere  are  built  into  the  model,  or  alternate  tables  of  atmosphere 
data  may  be  user  supplied.  Line  shape  options  include  Lorentzian,  gen- 
eralized Voigt,  and  collision  narrowed  line  profiles.  The  output  pro- 
vided by  this  model  are  the  absorption  coefficient,  Rayleigh  scattering 
coefficient,  and  transmittance  for  the  frequency  range  specified. 

Calculations  performed  by  the  attenuation  model  [4]  for  aerosols 
are  based  on  standard  Mie  theory  for  homogeneous  spheres  and  prov.,de 
the  scattering  and  absorption  coefficients  as  well  as  the  phase  function 
data  necessary  for  radiative  transfer  calculations.  Input  data  required 
are  the  laser  wavelength,  index-of-refraction,  size  distribution,  and 
particle  density  of  the  scattering  medium.  Several  size  distributions 
are  provided  as  "built-in"  options,  or  the  user  may  supply  slzt'  distri- 
bution and  particle  density  from  measured  data.  Output  from  this  model 
may  be  obtained  in  punched  card  form  for  use  as  input  into  radiative 
transfer  models  or  this  model  may  be  used  simply  as  a subroutine  In  the 
RT  model  as  in  the  present  application. 


IV.  RADIATION  TRANSPORT  MODELS 

To  predict  the  effects  of  radiation  absorption  and  scattering 
on  a laser  guidance  system,  a method  to  determine  the  amount  of  energy 
reaching  the  target  and  receiver  is  required.  Under  conditions  of  only 
moderate  aerosol  densities,  methods  which  consider  first-order  or  single 
scattering  effects  are  usually  adequate.  Under  conditions  of  very  dense 
aerosols  such  as  heavy  fog  or  smoke  in  the  transmitting  region,  however, 
second  and  higher  order  scattering  effects  become  important, and  more 
complex  procedures  are  required.  Models  using  both  methods  of  approach 
to  radiation  transfer  were  developed  during  this  effort.  An  existing 
Monte  Carlo  method  for  performing  multiple  scattering  calculations 
modified  and  adapted  to  the  present  problem  is  described  in  a separate 
report  [6].  The  remainder  of  this  section  presents  a method  employing 
first-order  scattering. 

The  geometry  for  the  model  presented  here  consists  of  a laser 
transmitter,  an  absorbing  and  scattering  medium,  the  target,  and  a 
receiver  viewing  the  target  (Figure  2).  The  target  receives  energy  by 
direct  transmission  and  by  scattering  from  aerosol  particles  located  in 
the  beam  and  reflects  part  of  the  energy  received.  Energy  reaches  the 
receiver  by  reflection  from  the  target  and  by  scattering  from  aerosol 
particles  in  the  transmitter  beam  located  within  its  field-of-view. 

The  transmitting  range  is  divided  into  segments  and  the  procedure  is 
used  to  calculate  the  energy  scattered  from  each  segment  to  the  target 
and  to  the  receiver.  The  target  area  is  divided  into  annular  rings  to 
determine  the  distribution  of  scattered  and  directly  transmitted  energy. 


i 
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RECEIVER/SEEKER 


The  basic  relationship  used  in  making  the  calculations  of  energy 
scattered  to  each  area  receiving  radiation  from  a segment  of  length  AZ 
becomes 


where 

* the  transmitter  output  power 

a = the  total  attenuation  coefficient  due  to  scattering  plus 
absorption 

" the  range  to  the  center  of  the  scattering  segment 

(jj  = the  albedo  or  ratio  of  scattering  to  total  extinction 

P(0j^)  ■ the  normalized  phase  function  at  the  angle  0^,  to  the 
receiving  area 

r^  ■>  the  path  distance  from  the  center  of  the  scattering  segment 
to  the  receiving  area  and  the  last  term  accounts  for 
attenuation  along  r^. 

The  total  energy  reaching  the  target  area  from  N segments  plus  directly 
transmitted  radiation  is 


P 


t 


P 

o 


N 


i-1 


The  scattered  energy  is  summed  over  all  the  annular  areas  of  the  target 
and  the  sum  over  N segments  gives  the  total  energy  scattered  to  the 
target  (neglecting  higher  order  scattering).  For  the  directly  trans- 
mitted radiation,  a gausslan  distribution  is  assumed  and  the  average 
intensity  is  computed  for  each  annular  area  from  scattered  and  directly 
transmitted  radiation. 

Reflected  energy  from  the  target  (assumed  to  be  a diffusely 
reflecting  surface)  that  reaches  the  receiver  is  computed  by  the 
relation 
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where 

p = the  target  reflectivity 

Cp  •*  the  receiver  viewing  angle  to  the  target  normal 

A = the  receiver  area 
r 

= the  range  from  the  target  to  receiver. 

For  a receiver  with  field-of-view  6^  viewing  the  target  at  angle  cp  to 
the  incident  beam,  the  beam  length  seen  by  the  receiver  is 


sin  7 


where  7«=it-<li-0/2  and  the  distance  from  the  receiver  to  the  beam 
' r 

intersection  point  Z is  given  by 

Li 

Z sin  CD 

Z - -E - 

L sin  7 

The  path  length  Z^  is  divided  into  M segments  of  length  the 

relation  used  to  determine  the  energy  reaching  the  receiver  from  all 
segments  becomes 


\hj 


where  7^  is  the  angle  (counterclockwise)  from  each  segment  to  the 

receiver  and  Z,  is  the  distance  from  the  center  of  each  segment  to  the 
Li 

receiver.  The  total  energy  reaching  the  recei’vr  by  reflection  and 
scattering  becomes 

P - P ^ + P 
r r,t  r,s 

In  practical  cases  of  interest  to  the  Army,  the  laser  transmitting 
region  consists  of  a clear  atmosphere  over  most  of  the  path  and  a smoke 
cloud  or  heavy  aerosol  concentration  over  the  remaining  distance.  Pro- 
vision has  been  made  in  the  model  described  here  by  the  introduction  of 
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a parameter,  giving  the  distance  from  the  transmitter  to  the 

beginning  of  the  aerosol  cloud.  Then,  in  performing  the  computations, 
if  no  absorption  or  scattering  takes  place.  When 

^ ^cloud’  computations  are  performed  as  previously  shown. 

Similarly,  in  performing  the  computation  of  energy  reflected  and  trans- 
mitted to  the  receiver,  no  attenuation  takes  place  once  the  cloud  has 
been  exited. 
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Appendix.  UTILIZATION  INSTRUCTIONS 


I 


I 


All  of  the  models  discussed  in  this  report  have  been  adapted  for 
use  on  the  Army  Missile  Materiel  Readiness  Command  (MIRCOM)  CDC  6600 
computer  at  Redstone  Arsenal.  This  appendix  describes  briefly  the  input 
data  required  for  each  of  the  computer  models  and  contains  Fortran 
listings  of  the  in-house  developed  models.  Instructions  for  use  of  the 
Monte  Carlo  radiation  transport  model  is  presented  in  a separate  report 
prepared  by  Radiation  Research  Associates  personnel  [7]. 

1.  Turbulence -Induced  Beam  Wander 


Required  input  data  for  the  turbulence- induced  beam  wander  model 
consists  of  run  options,  laser  transmitter  characteristics,  range  and 

2 

time  period,  tables  of  refractive  index  structure  constant,  C^^,  and 


windspeed  for  transmitter  and  seeker  paths,  array  sizes  for  frequency, 
and  range  segments.  The  following  list  contains  input  data  card  formats, 
variables  read  on  each  card,  and  a description  of  each  variable.  The 
CDC  6600  computer  program  library  must  be  attached  for  each  run. 


Card  Description 

1 lOPT,  NRUNS 
Format  (214) 

lOFT  >■  1,  calculations  for  designator  path  only 

■ 2,  calculations  for  designator  and  receiver  paths 
NRUNS  ■ Number  of  cases  to  be  run.  A complete  set  of 
the  remaining  data  cards  are  required  for  each 
case. 

2 LAMB,  DIAM,  THET,  TDOT,  RANG,  TIME 
Format  (6E10.4) 

LAMB  laser  wavelength  in  meters 

DIAM  diameter  of  laser  designator  aperture  in  meters 
THET  * beamspread  angle  of  laser  designator  in  radians 
TDOT  ■■  angular  slue  rate  of  laser  designator  in 
radians /second 

RANG  ~ distance  from  designator  to  target  in  meters 
TIME  *■  time  duration  of  beam  wander  calculation  in 
seconds 

3 N,  M,  N2 
Format  (3I4) 

N ■■  Number  of  segments  in  designator  path 
M " number  of  frequencies  in  power  spectrum  calculations 
N2  Number  of  segments  in  seeker  to  target  path.  Not 
needed  if  lOFT  •■1. 


ii 


y 

ii 


k 
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Card 

4 


Description 


CN(I),  I » 1,  N 

Format  (7E10.4)  2 _2/3 

CN(I)  "values  of  for  laser  designator  path  In  (meter)  ' , 

One  value  required  for  each  segment. 

5 V1(I),  I - 1.  N 
Format  (7E10.4) 

VI (I)  = values  of  crosswind  velocity  for  designator  path 
In  meters /second.  One  value  required  for  each 
segment. 

Note:  If  lOPT  = 1,  no  more  cards  are  needed.  If  10PT  = 2, 
the  following  cards  are  required. 

6 DIV,  RIV 
Format  (2E10,4) 

DIV  = seeker  aperture  diameter  In  meters 
RlV  " seeker  range  to  target  In  meters 

7 CN2(I),  I «=  1,  M 
Format  (7E10.4) 

CN2(I)  ■=  values  of  for  seeker  path  In  (meter)”^'^. 

One  value  required  for  each  segment. 

8 V2(I),  I " 1,  M 
Format  (7E10.4) 

V2(I)  " values  of  crosswind  velocity  In  seeker  path  In 
-2/3 

(meter)  . One  value  required  for  each 
segment. 


B 


2.  Molecular  Attenuation  Model 

The  following  list  presents  the  data  Input  needed  for  the  molecular 
absorption  model  In  use  at  the  US  Army  Missile  Research  and  Development 
Command  (MIRADCOM).  This  Information  was  extracted  from  the  report  by 
Miller  et  al.  [3j  which  Included  additional  Input  Information.  Spectro- 
scopic data  required  are  obtained  from  the  AFGL  tape  which  must  be 
requested  from  the  MIRCOM  tape  library  for  each  run.  The  current  serial 
number  of  the  AFGL  data  tape  used  Is  S08997. 

Card  Description 

1 ALTI,  ELVTN,  DIST,  IHR,  ICM,  IRL,  lAL,  ID 

Format  (3F5.1,  4ll,  6QX,  II) 

ALTI  " Initial  altitude  of  path  in  kilometers 
ELVTN  " angle  of  elevation  In  degrees  of  path  from  the 
horizontal 
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A 


r ^ 


Card 


L 


1 


2 


3 


Description 

DIST  = distance  of  path  in  kilometers 

IHR  * digit  1 to  include  high  resolution  molecular  line 
absorption  effects,  0 otherwise 
ICM  “ digit  1 to  include  molecular  continuum  absorption 
effects,  0 otherwise 

IRL  *=  digit  1 to  include  Rayleigh  scattering  effects, 

0 otherwise 

lAL  * not  used;  leave  blank  or  0 
ID  ■ digit  1 in  column  80 

VI,  V2,  CINC,  FINC,  SETBAK,  BOUND,  ACY,  ID 
Format  (2F12.4,  2F10.4,  2F8.4,  E10.3,  9X,  II) 

VI  “ starting  frequency  in  cm  ^ 

V2  * end  frequency  in  cm 

CINC  ■ coarse  increment  in  cm  ^ 

FINC  ■ fine  increment  in  cm  ^ 

SETBAK  ■ wave  number  interval  from  line  centers  at  which 
reversion  to  FINC  should  occur 
BOUND  ■ MAXIMUM  half  width  of  integration  range 
ACY  ■*  a measure  of  accuracy  desired  in  terms  of  magnitude 
ot  transmittance  (not  a 7.) 

ID  <■  digit:  2 (in  column  80) 

NGDN,  NADN,  ICPN,  IPRO,  NP,  DELV,  SLIT,  ID 
Format  (511,  F8,5,  F8.4,  58X,  II) 

NGDN  ■ number  of  atmospheric  model  or  gas  density  versus 
altitude  distribution 

0 (zero)  for  US  1962  Standard  Atmosphere 

■ 9 for  user  supplied  atmospheric  model 
NADN  * not  used;  leave  blank  or  0 

ICPN  • 0 for  continuum  calculation  only  at  one  wave- 
number  (the  average  value  of  VU  and  VL  will  be 
used  ordinarily) 

■ 1 for  continuum  calculation  at  every  high- 
resolution  wavenumber 

IPRO  ■ 0 for  Lorentzian  line  profiles 

■ 1 for  collisionally-narrowed  line  profiles 

■ 2 for  generalized  Voight  profile 
NP  ■■  not  used;  leave  blank  or  0 

DELV  ■ Increment  in  cm”^  for  frequencies  at  which  trans- 
mission results  are  degraded 
SLIT  ■ half  width  of  response  function  (triangular) 

SLIT  ■ 0,0  prevents  convolution  by  the  triangular 
slit 

ID  ■ digit;  3 (in  column  80) 
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Card  Description 

ALT(I),  P(I),  TEMP(I),  AIR(I),H20(I),  03(1)  , G8(I)  , G9(I)  , GLO(I)  , 
ID  Format  (F5.1,  E10.4,  F5.1,  E10.4,  E9.3,  E9.3,  3E10.4,  IX,  II) 

ALT(I)  ■=  altitude  in  kilometers  for  ith  data  set  for  gases 
P(I)  = pressure  in  mb  at  ALT (I) 

TEMP(I)  ■ temperature  in  degrees  Kelvin  at  ALT (I) 

AIR(I)  « density  of  air  in  g/m^  at  ALT(I) 

3 

H20(I)  = density  of  water  vapor  in  g/m  at  ALT(I) 

3 

03(1)  = density  of  ozone  in  g/m  at  ALT (I) 

G8(I),  G9(I) , GIO(I)  = densities  of  additional  optional 

3 

gases  in  mol/cm  at  ALT(I) 

ID  - digit  4 (column  80) 

5 DALT(I),  DENS (I),  ID 

Format  (F5.1,  E10.3,  64X  II) 

DALT(I)  = altitude  in  kilometers  for  ith  data  set  for 
aerosols  . 

DENS (I)  " density  of  aerosols  in  particles/cm  at  DALT(I) 
ID  “ digit  5 (column  80) 

A final  type  5 carrying  DALT(  ) = 999,0  is  required 
Note:  Type  4 and  5 are  needed  only  if  the  user  wishes  to 
insert  his  own  density  versus  altitude  models  for 
(Type  4)  gases  and/or  (Type  5)  aerosols.  Each  type 
requires  one  card  per  altitude. 

3.  Aerosol  Attenuation  Model 

Input  data  required  for  the  aerosol  attenuation  model  have  changed 
somewhat  from  the  format  presented  by  Gomez  et  al.  [4],  but  the  calcula- 
tions performed  remain  essentially  the  same.  Additional  size  distribu- 
tion options  have  been  built  into  the  model  and  the  input  data  have  been 
rearranged  slightly.  This  model  is  used  to  provide  the  phase  functions 
and  aerosol  attenuation  coefficients  for  the  First-Order  Radiation 
Transport  (FORT)  model. 


Card  Description 

1 WAVE , DENS 

Format  (2E12.6) 

WAVE  ■ laser  wavelength  in  micrometers 
DENS  ■ particle  density  in  particles /cc 


4 

(One  per 

j altitude) 


i 
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Card 


Description 


2 


f 


3 


IDSTP,  NRADI,  NCRDS,  IT,  MQRTE,  MCRTE 
Format  (615) 

IDSTP  “0,  1,  7,  size  distribution  option.  See 

Card  Type  3 for  details. 

NRADI  ■ number  of  different  radii  In  size  distribution. 
If  IDSTP  “ 5,  NRADI  is  used  to  give  number  of 
different  ranges  In  size  distribution. 

NCRDS  ■ 0,  for  printed  output  only 

* 1,  for  printed  and  punched  output 
IT  ■ number  of  terms  desired  In  the  phase  function 
expansion. 

MQRTE,  MCRTE  normally  blank,  but  can  be  used  to  obtain 
additional  printed  output. 

Data  to  be  Input  here  depend  on  the  value  of  IDSTP  given  on 
Card  Type  2. 


f 


t 

i 

I 


If  IDSTP"0  • "Arbitrary”  distribution 
F(I),  R(I),  I-l,  NRADI 

Format  (2E20.10),  number  of  Input  cards  required  ■ NRADI 
F(I)  > number  of  particles  of  radius  R(I) 

If  IDSTP"!  » log  normal  distribution 
RBAR,  SIGMA 
Format  (2E20.10) 

RBAR  " mean  radius 
SIGMA  " standard  deviation 
If  IDSTP"2  * Wynn/Dawes  exponential  distribution 
RLO,  RHI,  CUE,  A,  B 
Format  (5E12.6) 

RLO  ■ lower  radius  of  distribution  in  micrometers 
RHI  " upper  radius  in  micrometers 
CUE,  A,  B are  constants  in  distribution 
If  IDSTP-3  — Deirmendjian  Model  C. 

No  Input  needed;  parameters  are  fixed. 

If  IDSTP"4  — Junge  distribution 
RLO,  RHI,  CUE,  A 
Format  (4E10.4) 

RLO  " lower  radius  in  micrometers 
RHI  " upper  radius  in  micrometers 
CUE,  A are  constants  In  distribution 
If  IDSTP"5  • Modified  gamma  distribution 

RLO,  RHI,  RC,  ALF,  GAM,  DENS,  NRADI 
Format  (^12.6,  13) 

Number  of  cards  required  equals  value  of  NRADI  read  in 
Card  Type  2. 

RLO  " lower  radius  in  micrometers 
RHI  " upper  radius  in  micrometers 
RC  " mode  radius  in  micrometers 
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Card  Description 

3 ALF,  GAM  are  constants  in  distribution 
DENS  = particle  density  in  this  range 
NRADI  = number  of  radii  in  this  range 

i If  1DSTP=6  FOG  model 

VIS 

Format  (E20.10) 

! VIS  = visibility  in  kilometers 

' If  IDSTP=7  Hoidale  dust  model 

VIS 

Format  (E12.6) 

VIS  = visibility  in  kilometers 

4 EM,  CAY,  EMM,  CONG 
Format  (4F10.6) 

EM  = real  part  of  particle  index-of- refraction 
CAY  = magnitude  of  imaginary  part  of  particle  index-of- 
refraction 

EMM  •=  index-of-refraction  of  atmosphere  usually  1,0 
CONC  is  not  used  in  this  version:  internally  set  to  1.0. 

4.  First-Order  Radiation  Transport  (FORT)  Model 

The  aerosol  attenuation  model  described  in  the  previous  section  is 
contained  as  subroutines  in  the  FORT  model.  Data  input  for  the  FORT 
model  consist  of  the  three  cards  described  in  the  following  list  of  data 
followed  by  the  input  described  in  the  previous  section  for  the  aerosol 
model.  A complete  Fortran  listing  is  included  at  the  end  of  this 

Description 

LAMB,  DIAM,  THET,  RANG,  POWR,  ZCLOUD 
Format  (7E10.4) 

LAMB  ■ laser  wavelength  in  meters 
DIAM  * transmitter  aperture  diameter  in  meters 
THET  ■■  laser  beamspread  angle  in  radians 
RANG  ■ distance  from  transmitter  to  target  in  meters 
POWR  “ transmitter  output  power  in  meters 
ZCLOUD  ■■  distance  from  transmitter  to  aerosol  cloud  in 
meters 

ZR,  PHI,  THETR,  DR,  RHO 
Format  (5E10.4) 

ZR  ■ range  from  target  to  receiver  in  meters 
PHI  B viewing  angle  of  receiver  in  radians 
THETR  ■ field  of  view  of  receiver  in  radians 
DR  B diameter  of  receiver  aperture  in  meters 
RHO  B reflectivity  of  target 
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section. 

Card 

1 


Description 

3 NZ,  NR,  NM 

Format  (314) 

NZ  = number  of  segments  in  transmitter  path 
NR  = number  of  annular  areas  on  the  target 
NM  •=  number  of  segments  in  transmitter  path  within 
f ield-of-view  of  receiver. 

Card  Type  4,  etc.,  for  this  model  begin  with  Card  Type  1 described  in 
previous  section  of  this  appendix. 


FORTRAN  LISTING  OF  TlIRBULtNCE  INnuCED  BEAM  WANDER  MODEL 


■i 

\ 


PROGRAM  MAIN!  I NPUT  = 6'j  , OUT  PlJT  = 65,  TAPES  =I  NPIIT  , TAPE6  = OOT  PUT  ) 

C MlCOh  POINTING  JITTER  PROGRAM 

C calculates  riXBULENCE  INDUCED  POINTING  JITTER  AND  POWER  SPECTRUM 
C FOR  LASER  TARGET  DESIGNATOR  AND  TERMINAL  HOMING  SEEKER 

C »****RE0U1REU  INPUT  DATA»**»# 

C lOPT  = 1,  CALCULATIONS  FOR  DESIGNATOR  PATH  ONLY 

C 2,  CALCULATIONS  FOR  DESIGNATOR  AND  SEEKER  PATHS 

C NRUNS  = NO.  OF  CASES  TO  BE  CALCULATED  ( SEPARATE  SET  OF  DATA  IS 
C REUUIREO  FOR  EACH  CASEI  . 

C lamb  = LASER  WAVELENGTH  IN  MtTERS 

C DIAM  = LASER  TARGET  UESIGNATUR  APERTURE  DIAMETER  IN  METERS 

C THET  = LASER  BEAMSPREAD  ANGLE  IN  RADIANS 

C TOOT  = laser  beam  SLUE  RATE  IN  K AO | ANS/ SEC ONO 

C RANG  = PROPAGATION  RANGE  FROM  TARGET  DESIGNATOR  TO  SPOT  IN  METERS 

C time  = DURATION  OF  CALCULATIUN  OR  TEST  IN  SECONDS 

C CN(II  = values  of  REFRACTIVE  INOEX  STRUCTURE  CONSTANT  (CN)=i-*2 
C WITH  ONE  VALUE  EOR  EACH  SEGMENT  OE  RANGE  FROM  LASER 

C DESIGNATOR  TO  TARGET  (IN  ME T ERS** ( -2/ 3 ) I 

C Vl(I)  = SET  OF  VALUES  OF  CROSSWIND  VELOCITY  CORRESPONDING  TO 
C EACH  SEGMENT  OF  RANGE  FROM  LASER  DESIGNATOR  TO  T ARGE T I M / SE C ) 

CM  = NO.  OF  REOUENCieS  FlJR  WHICH  POINTING  JITTER  POWER  SPECTRUM 

C IS  TO  BE  CALCULATED 

C N = NO.  OF  SEGMENTS  OF  LENGTH  DELZ  FROM  DESIGNATOR  TO  TARGET 

C OlV  = DIAMETER  OF  SEEKER  APERTURE  in  METERS 

C RIV  = range  from  TARGET  TO  SEEKER  IN  METERS 

C CN2(I)  = values  of  REERACTIVE  index  STRUCTURE  CONSTANT  FOR  EACH 
C SEGMENT  OF  RANGE  FROM  TARGET  TO  SEEKER  ( MET ERS*» ( -2/ 3 ) ) 

C V2(l)  = VALUES  OF  CROSSWIND  VELOCITY  FOR  EACH  SEGMENT  OF  RANGE 
C FROM  TARGET  TO  SEEKER 

C N2  = NO.  OF  SEGMENTS  OF  LENGTH  OELIV  FROM  TARGET  TO  SEEKER  , 

dimension  CN(201,V1{20),FO(2U),RO(201,ER(1025),PS( 1025),V2(20) 
COMPLEX  RAN(20A8) 

COMMON/Z/  FO.RO.LAMB 

COMMON  /XX/  CN2(20),01V,R1V,N2 

COMMON  /YY/  RAN 

COMMON  til/  ORO 

REAL  LAMB 

EXTERNAL  DESUBfFALPH 
10  FORMAT ( 7E10. 4) 

20  FORMATOI^) 

READ(S,20)  IOPT, NRUNS  ' 

DO  1000  LU  = It NRUNS 

READ! b, 10)  LAMB. 01  AM, THET, TOOT, RANG, TIME 
REA0(b,20)  N,M,N2 
REAOIb.lO)  (CN(t),I  = 1,NI 
READ(b,10)  (V1(I),I  = 1,N) 

IF( IOPT  .EO.  1 ) GO  TO  1A5 
READIb.lO)  DIV.RIV 
READ(b,10l  (CN2( I ),I=1,N2I 


READ(t),10)  (V2(ri,£  = 1,N2I 

25  FORMAr(*  CALCULATION  OF  ROWtR  SRfcCIRdM  ANO  ruKbULENCE  INDUCED  FO I 
INTING  JITTER  OF  A LASER  TA-tGtT  DFS  I GNA  TOR  * // ) 

30  FORMAT!*  CALCULATION  OF  ROWER  SPECTRUM  AND  TURBULENCE  INDUCED  PO I 
INTING  JITTER  OF  A LASER  TARGET  DESIGNATOR  AND  SEEKER*//) 

^0  FORMAT!*  LASER  WAVELENGTH  =*,E1D.4,*  METERS,  DESIG.  APERT.  DIAM=* 
1,F10.6,*  METERS,  EEAMSPREAD  ANGLE  =*,F10.6,*  RADIANS*/) 

45  FORMAT!*  SEEKER  APERT.  DIAM.  =* , F 1 0.  f>,  *ME  T E RS  , RANGE  FROM  TARGET  T 
10  SEEKER  =*, F10.2,*METERS*/) 

60  FORMAT!*  BEAM  SLUE  RATE  =*  , E 10 .4 , *RAO/ SEC , DESIGNATION  RANGE  =*,F 
110.2,*  METERS*/) 

55  FORMAT!*  DURATION  OF  TEST  IS  *,F10.4,*  SECONDS*/) 

60  FORMAT!*  NO.  OF  SEGMENTS  IN  DESIGNATOR  PATH  =*,13/) 

65  FORMAT!*  NO.  OF  SEGMENTS  IN  SEE<tR  PATH  = *,13/) 

70  FORMAT!*  NO.  OF  FREQUENCIES  PUR  WHICH  POWER  SPECTRUM  IS  TO  BE  CAL 
ICULAT  EU  =*, I 4//) 

BO  FORMAT!*  VALUES  OF  RcFRACTlVF  INDEX  STRUCTURE  CONSTANT  AND  WIND  S 
IPEEU  IN  designator  PATH*/) 

85  FORMAT!*  VALUES  OF  REFRACTIVE  INDEX  STRUCTURE  CONSTANT  AND  WIND  S 
IPEED  IN  seeker  path*/) 

90  FORMAT!*  SEGMENT  NO.  = *,I2,*,  REF.  INDEX  STRUCTURE  CONST.  =*,£12 
1.6,*  !METER)2/3  , WIND  SPEED  =*,F10.4,*  METER/SEC*/) 

100  FORMAT!*  VALUES  OF  FREQUENCY  FOR  WHICH  POWER  SPECTRUM  CALCULATIO 
INS  ARE  TO  BE  MADE*/) 

1 10  format  I 6F16. 6) 

115  format ! 10E12. 4) 

120  FORMAT!*  FOR  FREQUENCY  OF  *,F10.4,*  HERTZ,  THE  CALCULATED  POWER  S 
IPECTRUM  IS  *E12.4,*,  THE  VARIANCE  1S*E12,4/) 

130  FORMAT! IHO) 
l4o  format !1H1) 

145  CONTINUE 

WRIT  E!  6,  140  ) 

IF!IOPT  .EO.  1)  WRITE!  6,25) 

IF!IOPT  .£0.  2)  WRITE!  6,30) 

WR1TE!6,40)  LAMB, DIAM, THET 
IFIIOPT  .EO.  2)  WRITE!  6,45)  D1V,R1V 
WRITE!6,50)  TDOT.RANG 
WRITEI6,55)  TIME 
WRITE!6,60)  N 

IF!IOPT  .EO.  2)  WKITEI6,65)  N2 
WRITE!6,70)  M 
WRITE! 6,80) 

DO  150  I = 1,N 

150  WRITE! 6,90)  I , CN ! I ) , V 1 ! I ) 

WRir£|6,130) 

IF! lOPT  .Eg.  1)  GO  TO  156 
WRITE!6,85) 

DO  155  I = 1,N2 

155  WRITE!6,90)  I , C N2 ! I ) , V2 1 I ) 

WRITE! 6, 130) 

COMPUTATION  OF  TIME,  FREQUENCY  AND  SPATIAL  INCREMENTS 

156  CONTINUE 

DELT  = TIME/M  »DELF  = l./TIME 

DELZ  = RANG/FLOAT! N)  IDFLIV  = < 1 V/ FLOA T! N2  ) 
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MM  = M + M $M1  = M + 1 SMMl  = MM  + 1 1MM2  = MM  + ? 

MSO  = SORT( FLQAT(MM) 1 
DO  160  1 = 2«M1 
FR<I)  = (t  - 1)*0ELF 
160  CONT  INUE 

R2  = OIAM/THET 

R - RANG  * R2 

02  = OIAM  + THET*RANG 

C COMPUrAriON  of  effective  wind  velocity,  coherence  length  ANO 
c normalization  frequency  for  each  segment  of  path  from  laser 
c designator  to  target 

21  = R2 
Z1  = UEL2/2. 

ROT  = 0.0 
DEL  = 0.0 
DO  200  I = 1,N 
21  = ZI  + Z1 
Z1  = UEL2 

VEI  = Vl( 1)  ♦ TOOT*(ZI  - K2) 

ROin  = ( 16.7*DELZ*CN(  n*IZl/R  )♦*  1 . 6666  7/ ( LAMB*LAMB  ) ) 

FOin  = VEI  /(3.16159>»D2*ZI/R) 

BCT  = ROT  + R0( I ) 

DEL  = UEL  + OELZ*CN(  n*(RANG  - ZD/RANG 
R0( I I = R0( I )**( -.6) 

200  CONTINUE 

ROT  a K0T**(-.6) 

C COMPUTATION  OF  BEAM  SPREAD  ANGLE  DUE  TO  TURBULENCE  ANO  OIFFRACTION 
C computation  of  spot  diameter  ON  TARGET 
DRO  « DIAM/ROT 

CALL  IGRATIO., 1. , .0 1 , 1 , OESUB,RDRO ) 

RDRO  = 1.0/( SORT! 5.092958*(OKO)**2*RORO ) ) 

THET12  = 1. 128*LAMB/ROT*RORO 
OTHET  = S0RT(THET12**2  + THET**2) 

022=  UIAM  + DTHET^RANG 
WRirE(6,210)  D2,R2,R,OEL,D22 
WRITE! 6,240) 

WRITE(6,110)  (R0(I),I=1,N) 

WRITE(6, 130) 

WRirE(6,245) 

WRITE(6,1I0)  (F0(I),I=1,N) 

210  FORMAT!*  SPOT  OIAM.  =*,E12.6,*,  R2  =*,F10.A,*,  R1+R2  =*,F10.4,*, 
1 DEL  = *,E12.6/,*  TURB.  INDUCED  SPOT  OIAM.  = *E12.6/) 

240  FORMAT!*  VALUES  OF  RO!I),  I = 1,  N)*/) 

245  FORMAT!*  VALUES  OF  FO!I),.I  = 1,  N)*/) 

C COMPUTATION  OF  ANGLE  OF  ARRIVAL  POWER  SPECTRUM  OF  LASER  DESIGNATOR 
F2  = 0. 

PS! 1 ) = 0. 

DO  300  J = 2, Ml 
F = FR!J) 

FI  =0. 

CALL  SPECT!F,F1,02,N) 

PS!J)  = FI 

IF!IOPT  .EO.  1)PS!J)  = F1*!D2/DIAM)**2 
F2  = F2  ♦ PS!J)*OELF 
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IF(J  .LT.  M»GQ  TO  300 

IFIIOFT  .EO.  nwRITE(6,  120IF,PS(J)tF2 
300  CONTINUE 

IF( lOPT  .EO.  n GO  TO  306 

C computation  of  effective  wind  VELOCITY,  COHERENCE  LENGTH  AND 
C normalization  FREOIIENCY  FOR  EACH  SEGMENT  DF  PATH  FROM 
C TARGET  TO  SEE<ER 

21  = 0.  $21  = 06L1V/2. 

DO  220  I = 1,N2 
21  = 21  + 21 
21  = UELlV 

VEI  = V2(I)  + TDOT*(RANG  - 21  I 

R0( 1 ) = ( 16.700  ELlVOCN2( 1 )o( 2 I /K 1 V ) 00 1 . 6666 7/ t LAMB0LAM3 ) )«0( -.6 ) 
FOd)  = VEI/(3.l415qoniVo2I/KlV) 

220  CONT  INUE 

WR  ITE(  6,  130  ) 

WRITe(6,2400) 

WRITE!  6,  110)  (R()(  1 ) , I = 1,N2) 

WRir6(6,130) 

WRITE(6,2A50) 

WRITE(6,110)  (Foil ), 1=1, N2) 

2A00  FORMAT!*  VALUES  OF  R!)!I),  I = 1,  N2)*/) 

2^50  FORMAT!*  VALUES  OF  FO ! I ) , I = 1,  N2)*/) 

C computation  of  TURBULENCE  INDUCED  POINTING  JITTER  POWER  SPECTRUM 
C FROM  TARGET  SPOT  TO  LASER  SEEKER.  COMPUTATION  OF  TOTAL  POWER 
C SPECTRUM  FROM  LASER  DESIGNATOR  TO  SE6<ER  AND  POWER  SPECTRUM  VARIANCE 
CALL  THETOITHETAO.FALPH) 

WRITEI6, 140) 

F2  = 0. 

DO  305  J = 2, Ml 
F = FR!J) 

FI  = 0. 

CALL  SPECT I F, Fl,01V,N2 ) 

PSIJ)  = PS!J)  + FI/! 1. +! 02/ (RlV*rHETAO) ) **2 ) 

PSIJ)  = PSIJ)*I02/01V)**2 
F2  = F2  + PS(J)*DELF 
I F( J .LT.  M)G0  to  305 
WRITE!6,120)  F,PSIJ),F2 
305-  CONTINUE 
306  CONT  1 NUE 

WRITE!6,310) 

WRITE!6,115)  !PS!J),J  = 1,M1) 

WRITE!6,130) 

310  FORMAT!*  CALCULATED  Pf)WER  SPECTRUM  VS.  FREOIIENCY*/) 

DO  1000  L = 1,2 
WRITE!6,140) 

WRITE!6,320)  L 

320  FORMAT!*  /////  OUTPUT  Ff)K  DIRECriTjN  *,I2,*  /////*//) 

C generation  of  random  SEUUENCE  HAVING  SAME  POWER  SPECTRUM  VARIANCE 
C AS  INDUCED  BY  TURBULENCE.  ADD  SYMMETRIC  TERMS  FOR  NEGATIVE 
C FREOUENCIES.  COMPUTE  MEAN  AND  VARIANCE  OF  RANDOM  ARRAY 
RANIl)  = 10. ,0.) 

DO  350  I = 2, Ml 
MMM  = HM2  - I 


II 
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RAN(I)  - DNRHALlO.t 1.0)«SORT(PS(  I)/OELT) 

RAN(MMM)  = RAN( I ) 

3 50  CONTINUE 

WRITE! 6, 130) 

c compute  and  write  mean  and  variance  oe  ranooh  array 
WRITE(6,360) 

CALL  HEANVARI ItM) 

C FAST  FOURIER  TRANSFORM  RANDOM  ARRAY 
CALL  FFT(RAN,MM,+1 ) 

DO  ^50  I = 1,MM 
RANI  I ) = RANI  I )/MSO 
450  CONTINUE 

C COMPUTE  AND  WRITE  MEAN  AND  VARIANCE  OF  TIME  SEOHENCE. 

WRIT  El  6,480) 

CALL  MEANVARI M1,MM ) 

WRITEI6,130) 

C WRITE  TRANSFORMED  ARRAY  VALUES  CORRESPOND  INF,  TO  TIME  VALUES 
C OF  POINTING  JITTER  FOR  ONE  DIRECTION. 

WRITEI6,130) 

WRITEI6,460)  OELT 

WRITEI6,115)  IRFALIRANII )),I  =M1,MM) 

1000  CONTINUE 

360  FORMATI*  mean  ANO  VARIANCE  OF  RANDOM  ARRAY*//) 

460  FORMATI*  VALUES  OF  POINTING  JITTER  AT  *,F10.6,*  SEC  INTERVALS  BEG 
lINNING  AT  T = 0*/) 

470  FORMATI*  RANDOM  VALUES  VS.  FREOiJENCY  AT*,F10.4,*  HZ  INTERVALS*/) 
480  FORMATI*  MEAN  AND  VARIANCE  OF  TIME  SEQUENCE*//) 

END 

SUBROUTINE  SPECTI F, FI ,02,N) 

COMMON/Z/  F0I20),R0I20),LAMB 

real  lamb 

FACT  = 1-32E-2*! LAMB/D2)**2 
FI  =0. 

DO  250  I = 1,N 

IFIF  .LE..332*F0II ))  G =1. 

IFIF  .GT..332*F0I I ) ) G =1.12  - .361*F/F0II) 

. IFIF  .GE.3.10*F0I I ) ) G =0. 

FI  = FI  + FACT*I I 02/R0I I ) )**5/IF*F*F0I  I ) ) )** .33333*G 
2 50  CONTINUE 
RETURN 
END 

SU3R0UTINE  OESUB I X , Y , NOEO ) 

COMMON  /ZZ/  DRO 

1 V=X*I I ACOSIX)-X*l l.-X**2)**.5)*EXPI-3.44*I  DRO*X ) ** 1 . 666 7* I l.-X**0. 
13333) ) ) 

RETURN 

END 

SL8R0UTINE  TH ETO I THETAO , FAL PM) 

COMMON  /XX/  CN2I 20) ,01V,R1V,N2 
DELIV  = R1V/FL0ATIN2) 

C calculate  OlINF 

011NF=0.  $Sl=0ELlV/2.  $S  =0. 

D1V3  = IDlV)**I-.3333) 

DO  20  I = 1,N2 
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S = S + SI 

D1INF  = DEL1V*CN2(  n*(  ( S/«1V)**1.6^67)  * OlINF 
20  SI  = OELIV 

OlINF  = 0.5*11.97*DlV3*ni£NF 
C INITIAL  ESiIHATE  FOR  THETAO 
WRITE(6,35)  OlINF 
THETAO  = l.E-4 
2 5 CONTINUE 

XU=0.  $DirHE  = 0.  $S  = 0. 

SI  = UELlV/2. 

DO  30  I = 1,N2 
S = S + SI 

Xll  = THErAO»(RlV-S)/r)lV 

DITHE  = 0EL1V*CN2(I  )*(  (S/RlV)«<>1.6ftft7)*‘FALRH(XIll  + (jlTHF 
30  SI  = UELIV 

DITHE  = U1THE*D1V3 

I F(  ABS(  ( 01 INF-OITHE  ) /01  INF  ) .LT.  .001)  GO  Tn<*0 
THETAO  = THETA0«(1.  + . 5* ( 01  £ NF-0 1 THE ) /0 1 1 NF ) 

WftITE{6,35)  THETAOtOlTHE 
GO  TO  25 

35  F0RMAT(2E16.a ) 

AO  CONT  I NUE 

WRITE(6,35)  OITHE 

RETURN 

END 

FUNCTION  FALPHIXI  1) 

DIMENSION  A(8) 

DATA  ( A< I ) , Ial,a )/1.98  71A, 10.  3433 1 -5. 90301  » 1 . 8361 9 , -0 .30 1 AA2 * 
1 2.  51  509E-2, --9. 752 29E-4, 1.356186-5/ 

IFIXIl  .GE.  ,5623)  GO  TO  10 
FALRH  = 10.66*1 (XI 1)**2) 

GO  TO  40 

10  IFIXIl  -.07.31.62)  GO  TO  20 
FALPH  = 0.0 
DO  15  I = 1,8 

15  FALPH  = A(I I*((XI1)**(I-1)I  * FALPH 
GO  TO  40 

20  IFIXIl  .GT.IOOO. ) GO  TO  30 
FALPH  = 7.8*1 (X  I 1 )**.06  ) 

GO  TO  40 

30  FALPH  = 11.97 
40  RETURN 
END 

SUBROUTINE  MEANVAR I N1 , N2 ) 

COMPLEX  RANI  2048) 

common  /yv/  ran 

REAL  MEAN1,MEAN2 

MEAN1=MEAN2=VAR1=VAR2=0. 

Do  400  I = N1,N2 
MEANl  = KEALIRANII))  -f  MEANl 
MEAN2  =:  AIMAGIRANI  I ) ) * MEAN2 
VARl  = (REALIKANI I ) ) )**2  + VAHl 
VAR2  = (AIMAGIRANI I )) )**2  + VAR2 
400  CONTINUE 


27 





MEANl  * MeANl/FL0AT(N2  - Nil 
MEAN2  > MEAN2/FL0AT(  N2  - Nl) 

VAHl  » VAK1/FU0AT(N2  - Nl I 
VAR2  » VA»<2/FL0AT(N2  - Nil 
WRirE(6.700IMEANl,  MEAN2 
WR1TE(6,800I VARlt VAR2 
RETURN 

700  FORMAT!*  MEAN  OF  REAL  PART  =*,E12.6,*,  MEAN  OF  IMAG  PART  =*£12.6/) 

800  FORMAT!*  VAR.  OF  REAL  PART  =*,E12.6,*,  VAR.  OF  IMAG  PART  =*E12.6/) 
END 
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1 

PROGRAM  FORr(  INPUT  = 65,0ljrPl/T  = A5,  IAPe5=lNPUT,  T A Pe6=0U f PUT  ) 

j 

c 

PROGRAM  TO  COMPUTE  FIRST  OROEK  RAOIATION  TRANSPORT  FROM  A LASER 

000110 

1 

c 

transmitter  to  a target  and  the  reflected  and  SCATTERED  RADIATION 

000120 

t 

t 

c 

TO  A RECEIVER  VIEWING  THE  TARGET 

000130 

■ ' 

c 

♦♦♦♦♦reuuireu  input  oata****» 

000140 

1 

c 

lamb  = LASER  WAVELENGTH  IN  MICROMETERS 

000150 

1 

1 

c 

DIAM  = TRANSMITTER  APERTURE  DIAMETER  IN  METERS 

000140 

1 

c 

THET  = laser  BEAMSPREAD  ANGlE  IN  RADIANS 

000170 

■ 

c 

RANG  = PROPAGATION  RANGE  FROM  TRANSMITTER  TO  SPOT  IN  METERS 

OOOIRO 

c 

POWR  = TRANSMITTER  OUTPUT  POWER  IN  WATTS 

000190 

r 

ZCLOUD  = DISTANCF  FROM  TRANSMITTER  TO  CLOUD  IN  METERS 

000200 

c 

ZR  = RANGE  FROM  TARt'.ET  TO  RECEIVER  IN  METERS 

000210 

c 

PHI  = VIEWING  ANGLE  OF  RECEIVER  IN  RADIANS 

000220 

c 

THETR  = FIELD  OF  VIEW  OF  RECEIVER  IN  RADIANS 

000230 

>3 

DR  = receiver  diameter  in  meters 

000240 

c 

KHO  = REFLECTIVITY  OF  TARGET 

000250 

c 

NZ  = NO.  OF  segments  0-=  LENGTH  DELZ  FR/)M  TRANSMITTER  TO  TARGET 

0002  40 

; 

c 

NR  = NUMBER  OF  ANNULAR  AREAS  IN  THE  TARGET  SPOT 

000270 

■ ! 

c 

NM  = number  of  segments  of  TRANSMITTER  PATH  SEEN  BY  RECEIVER 

000280 

• c 

***•*♦****«:*«♦♦•*•****♦**•**»#***»*•♦♦******** ************#**»******O0D2Q0 

Dl MENS  I ON  ZN(20),RN( 20) , PI 10U),C I 100) .AREA  120) ,PWHI20 ) ,AINTI20  ) 

000300 

DIMENSION  GAM(20).ZL( 20) 

000310 

C0MM0N/B<2  / C.  ALBDO.LLLL.NCROS.IT.ITT.NRADI 

000320 

REAL  LAMB.KEXT 

000330 

10  FOHMAT(7E10.4) 

000340 

20  format (31  A) 

000350 

READ! 5, 10)  LAMB.OI AM, THET, RANG, POWR, ZCLOUD 

000360 

READ(5,10)  ZR, phi , THETR, DR, RHO 

000370 

«EAD(i>,20)  NZ,NR,NM 

000380 

AO  FORMAT!*  LASER  WAVELENGTH  =*,F10.A,*  MICRON,  OESIG.  APERT.  DIAM= 

*000390 

l,Fl().6,*  METERS,  BEAMSPREAD  ANGLE  =*,F10.B,*  RADIANS*/) 

000400 

, 

50  FORMAT!*  DESIGNATOR  RANGE  =*,F10.2,*  METERS*/,*  DESIGNATOR  POWEROOOAIO 

1 =*,F10.3,*  WATTS*/,*  DISTANCE  TO  CLOUD  =*,FI0.2,*  METERS*/) 

000420 

55  FORMAT!*  RECEIVER  RANGE  =*,F10.2,*  METERS,  REC.  FOV  =*,F10.5,* 

R000430 

IADIANS*/,*  receiver  viewing  angle  =*,F10.A,*  radians,  receiver  niOOOAAO 

2AM.  =*,F10.6,*  METERS*/,*  TARGET  REFLECTIVITY  =*,F10.6/) 

000450 

60  FORMAT!*  NO.  OF  SEGMENTS  IN  DESIGNATOR  PATH  =*,I3/*  NO.  OF  ANNUL000A60 

lAR  elements  in  SPOT  =*,I3/,*  NO.  OF  SCATTERING  SEGMENTS  =*,13/1 

000470 

70  F0RMATI8F16.8) 

000480 

115  format ( 10E12.A1 

000490 

130  FORMAT(IHO) 

000500 

140  format (IHl) 

000510 

PI=3. 141 592653 5898 

000520 

WRITE! 6, 140) 

000530 

WR1TE(6,40)  LAMB, DIAM, THET 

000540 

WRITE! 6,50)  RANG , POWR , Z CLOUD 

000550 

WRIT E! 6, 55) ZR, THETR, PHI , OK, RHO 

000540 

WRITE!6,60)  NZ,NK,NM 

000570 

c 

COMPUTATION  OF  SPATIAL  INCREMENTS 

000580 

c 

COMPUTATION  OF  SPOT  DIAMETER  ON  TARGET 

000590 

IT  = NR*NZ 

000  600 

I OUT  = IT  - 2 

000610 

ZD  = UIAM/THET 

000620 

■ 

■ 
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ZC  = ZO  ♦ ZCLOUD 

000IS30 

ZT  < RANG  * 20 

000  W,0 

OT  s OIAM  > THET«RANG 

OOOftSO 

ZI  » ZC 

OOOftiSO 

OELZ  » (ZT  - ZCl/FLOATINZ) 

OOOf.70 

DEU  « OT/FLOAT(NRI 

000630 

21  » OELZ/2. 

000690 

00  200  I < l.NZ 

000700 

21  « ZI  ♦ Zl 

000710 

Zl  * OELZ 

000720 

ZNd)  > ZI 

000730 

200 

CONTINUE 

0007A0 

WRITEI6.210)  OT.Zn.ZT 

0007S0 

2 10 

FORMAT!*  SPOT  f)IAM.  =*,E12.6,*t  ZD  =»,F10.A,*,  ZT  = *,F10.A/) 

000760 

WRITE(6,215) 

000770 

WRITE(6,70)  (ZNIITiI  » l.NZ) 

0007fl0 

215 

FORMAT!*  VALUES  OF  (ZNd), I = l.NZ)  */) 

000790 

c 

:oMpurArioN  of  radius  to  annular  elements  dn  target 

OOOROO 

c 

COMPUTATION  OF  AREA  OF  ANNULAH  ELEMENTS  ON  TARGET 

000810 

RI  > 0. 

000820 

AI  - 0. 

000830 

Rl  = OELR/2. 

000860 

00  220  I » 1,NR 

000850 

PWRd  ) * 0. 

000860 

Al  ■ Pr*dOELR*FLOATd  I 1**2) 

000870 

AREAd)  * Al  - AI 

000880 

AI  « Al 

000890 

RI  ■ Rl  * Rl 

000900 

Rl  - OELR 

000910 

RNd  ) ■ RI 

000920 

220 

CONTINUE 

000930 

MRirE(6,130> 

000960 

WRITE!6,225) 

000950 

225 

FORMAT!*  VALUES  OF  ! (RNd  ), AREAd  ))  , I = 1,NR)  •/) 

000960 

MRITE!6,70)  ! (RN!  I ) , AREAd  ) ) , I * 1,NR) 

000970 

00  230  J ■ l.NZ 

000980 

DO  230  1 - l.NR 

000990 

ANG  * ATAN!RN!( )/!ZT  - ZN!J))) 

001000 

L « NR*!J  - 1)  ♦! 

001010 

230 

C(L)  = COS(ANG) 

001020 

KKK  - 1 

001030 

CALL  AGAUSS!P,KEXT.KKK) 

WRITE(6,  130) 

001050 

c 

WRITE16,235) 

001060 

c 

235 

FORMAT!*  VALUES  OF  !!C!I),  P!I)),  I = l.ID  */ ) 

001070 

c 

MR  IT  El  6, 70)  (Cin.Pd  ).C!  1*1 1.P!  1*1  ),C!  1*2).  P!  1*2  ).C!  1-^3  ).P!  1*3), 

001080 

c 

1I-1.I0UT,4) 

001090 

00  2^0  J - l.NZ 

OOllOO 

ZTZN  ■ ZT  - ZN!J) 

OOlllO 

POMFAC  ■ POWR*EXP!-KEXT*!ZNIJ)-ZC) )*! l.-EXP!-KEXT*OELZ ) ) 

001120 

1*EXP(-KEX T*2T2N)*ALBOO 

001130 

00  2^  I - 1,NR 

001160 

L ■ NR*!J  - 1)  *I 

00 1 1 50 

PWRIJ)  - PMRIJ)  * PnNFAC«P!L>*ARFAd)/(ZTZN**2  *RN!I)**2) 

001160 

30 


2 40  CONTINUE 

WftITE(6, 130) 

WRITE! 6,245) 

245  format  (•  ROWER  ANO  AVERAGE 
TPOW  = 0.0 

POWER  = POWR*EXP(-KEXT*(RANG 
DO  250  I = 1,NR 
TPOW  = TPOW  + PWRI I ) 

A INTI  I ) = PWR  ( I ) / AREA!  I ) 

2 50  CONTINUE 

WRIT  El  6, 70)  ((PWR(I),AINT(1  )),I=1,NR) 
WRITEI6,130) 

WRITE! 6,255)  POWER, TPOW 
2 55  FORMAT!*  DIRECT  POW-iR  TO 
1 POWER  TO  target  =*F16.8,» 

DO  2 60  I = 1 ,NR 
AI  NT  ! I ) = A1  NT!  I ) -f  2. 

1)  ) 

2 60  CONTINUE 

WRITE!6,265) 

WRITE!6,70I  ! AINT! I ) , 1=1,NR ) 

265  FORMAT!*  TOTAL  INTENSITY  = DIRECT 
WRITE!6,130) 

THETR2  = rHETR/2. 

GAMA  = PI  - THETR2  - PHI 
AR  = PI*0R*DR/4. 

2M  = 2R*SIN!THETR2I/SIN!GAMA) 

IFI2M  .GT.  RANG)  ZM  = RANG 
ZS  = RANG  - ZM 

IFIZS  .LT.  ZCLOUO)  ZS  = ZCLOUO 
IT  = NM 
lOUT  = IT  - 3 
IF! lOUT  .LT.  1)  lOUT  = 1 
OELZ  = IRANG  - ZS) / FLOAT! NM) 
WRITE!6,267)  GAMA,ZM,ZS 
267  FORMAT!*  GAMA  =*,F10.5,*,  ZM 
ZI  = ZS 
Z1  = UELZ/2. 

DO  270  I = 1,NM 
ZI  = ZI  + Zl 
Zl  = UELZ 
ZN!I  ) = ZI 


INTENSITY  IN  ANNULAR  AREAS  */) 


ZCLOUO) ) 


TARGET  =*,F16.0,*  WATTS*/,* 
WATTS*/) 


001170 
00 11 80 
001190 
001200 
001210 
001220 
001230 
001240 
001250 
001260 
001270 
0012R0 
001290 

SCATreRED001300  . 

001310 
001320 

*PnwER/! PI »Dr*Dr )*EXP !-2. *RN I 1 )*RN ! [ ) / I 0 T*DTOO I 330 

001340 
001350 
001360 
001370 
0013R0 
001390 
001400 
001410 
001420 
001430 
001440 
001450 
001460 
001470 
001480 
001490 
001500 
001510 
001520 
001530 
001540 
001550 
001560 
001570 
0015R0 


+ SCATTERED  VS.  RN  ! I )*/  ) 


*,F10.3,*,  ZS  =*,F10.3/) 


ZL!I)  = S0RT!ZR**2  + ! R ANG-ZN! I ) )**2  - 2 .*ZR* ! RANG-ZN ! I ) ) *COS ! PHI ) )00 1 590 


TETAR  = ASINIIRANG  - ZN ! I I )* SIN! PHI ) /ZL ! I ) ) 

GAM! I ) = PI  - TETAR  - PHI 
C!I)  = COSIGAMID) 

270  CONTINUE 

WRITE! 6,275) 

275  FORMAT!*  VALUES  OF  ! ! ZN! I ) ,ZL ! I ) ) , I = 1,NM)*,/) 
WRITE! 6,70)  ! ! Z N! I ) , ZL ! I ) ) , 1 * 1 , NM) 

WRITE!6,130) 

KKK  = 2 

CALL  AGAUSS!P,KEXT,KKK) 

PWR SR  = 0.0 


001 600 
001610 
001620 
001630 
001640 
001650 
001660 
001670 
0016R0 

001700 
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00  2S0  J 3 1,NM  001710 

ZLJ  ■ ZL(J)  0017?0 

IFIZCLOUD  .GT.  .001  .AND.  ZK*COS< PHI!  .GT.  ( RANG-ZCLfJUO ) I 001730 

IZLJ  = CZNtJ)  - ZCLa)0)/C0S(PI  - GAM(J))  001740 

PWRSR  » fWRSH  >*0WR*EXP(-KexT*(2N<  J)-2C)  )*l  1 .-EXP(-KEXr*0Ei.2  ) )*  0017S0 

1ALB00*P(J)*AR/(ZL(  J)**2)«EXP(-KEXT«(ZI.J)  ) 00 17 AO 

280  CONTINUE  001770 

ZRCL  3 ZR  0017R0 

IFIZCLOUD  .GT.  .001  .ANO.  ZR^COSIPHl)  .GT.  ( RANG-ZCLODO ) ) 001790 

IZRCL  3 (RANG  - ZCL QUO) /COS (PHI)  OOlfiOO 

PWRTR  3 (POWER  ♦ TPOW)*RHO*COS(PHn*AR/(  ZR*ZR*P1  )*EXP(-i<EXT*ZRCL)  001810 
WRITE(6.28S)  PWRSKtPWKTK  001820 

WRITE! 6. 130)  001830 

POWER  SCATTEREO  PRUM  BEAM  TO  RECEIVER  =*,E1A.8,*  WAT T S* /OO 1840 


285  FORMAT!* 

1,*  POWER  REFLECTED  FROM  TARGET  TO  RECEIVER  =*,E16.H,^  WATTS*/) 

END 

SU3R0UTINE  AG AUSS ( P SUM, UTSUM, KKK) 

DIMENSION  F(  500), R(  500 ) , FSUMI ( 1 O) ,C ( 100 ) ,P ! 1 00 ) , PSUM ( 100 ) 

C0MM0N/BK2  / C,  ALBOO, LLLL ,NCROS . 1 T , I TT ,NR AD  1 

1 format (615, 49X,I1) 

2 format (4F10. 6) 

3 FORMAT(2^X,2(  E20. 10) ) 

4 format (2E20. 10) 

5 F0RHAT(6E12.6,13) 

6 format (4E10. 4) 

7 FORMAT! IH  , lOOH****** AEROSOL  DISTRIBUTION  TYPE  IS  UNDEF I NeD*******00 19A0 

*EXECUriON  CONTINUING  ASSUMING  NO  AEROSOL  MATERIAL)  001970 

8 FORMATdH  ,41H  DISTRIBUTION  WAVELENGTH  REFRACT  I VE 1 1 3X , 12HC  E001980 

♦XTINCTIONi 13X,12HC  SCATTER  I NGt 16X,5HALBOO/ IH  ,6X , AHT YPE ,6X , 9H (M I CROO 1990 
*ONS) t8X,5HINOEX,l6X,12H(SO  M I CRONS ) 1 1 3X , 12H ( SO  MICR0NS)/1H  ,I9,4X  002000 
*f  F11.4,F10.4,2HH,F7.4,2HI  ) ,3E25.  14) 

9 FORMAT! lH0f4( 5H  MUt 8 X, 17HP HASE  FUNCTION  )/lH  ) 

10  FORMATdH  , 4(  FI  3.9,  E17.  10 ) ) 

11  FORMAT! 1H0,65X, 14HPHASE  FUNCTION/IH  ) 

12  FORMATdH  , 5X  , IHL,  20X , 16HL-TH  COEFF IC  1 ENT  ,2  3X  , 14MRMS  DEV  I AT!  ON /1H002050 

002080 

,8H  MU  ,4!30H  ORIGINAL  EXPANDED 

,F8.5,8E15.8) 

,/////24H  AEROSOL  PARAMETERS  ARE  ) 
format !1H*,24X,6HRBAR3  ,E20.  10, 1 5X , 7HS IGMA3  ,E20.10/) 
format  !1H'»,24X,5HRL0«  , E10.4,  1X,5HRH13  , E 10.4 , 1 X , 5HCUE3 
«3HA-  ,E10.4,1X,3HB>  ,E10.4/) 

18  format  !lH'f,24X,5HRL0-  , E 10. 4,  1 X,5HRH  I 
•3HA3  ,E10.4, 1X,4HVIS-  ,E10.4/) 

19  format  !lH-f,24X,5HKL0-  , E 10. 4,  1X,5HRHI  3 , E 10 .4 , 1 X ,4HRC  > 


001850 

001880 

001880 

001890 

001900 

001910 

001920 

001930 

001940 

001950 


002010 

002020 

002030 

002040 


* ) 

13  FORMATdH 

14  FORMAT!  IH 

15  FORMAT  ! IH 

16 
17 


)/lH  ) 002070 

0020RO 
002090 
002100 
,E10.4,1X,002110 
002120 
,E10.4, IX, 002130 
002140 
,E  10.4,  IX, 500  2150 

«HALF3  ,E10.4, 1X,5HGAM3  ,E10.4,7H  NRAD|3,I3,4H  JZ3,I3/)  002180 

20  format !1H*,24X,3HA«  , F4. 1, 1 X, 7HNRAQI 3 . 14 , 1 X ,5MGAM3  ,F4 . 1 , IX ,5HV I S0021 70 

*3  ,E2U.10/)  002180 

21  format !1H*,24X,40H  RADIUS  (MICRONS)  RELATIVE  DENSITY  /)  002190 


,E10.4,1X,5HCUE3 


22  FORMAT!*  INOX-*, 13,*  M-  «,F1U.6,*  K » -*,F10.6,*I.  EMM3*,P8.6,* 
*CONCENTRATION  - «,f  10.6  ) 

23  FORMAT!*  THIS  IS  A MIXED  CASE SUBSEQUENT  REFRACTIVE  INDEX 

*OUTS  ARE  NOT  GENERALLY  VALID*) 

24  FORMATdH  ,24X,«  HIE  SIZE  PARAMETER  RANGE 


002200 
002210 
PKINT-002220 
002230 
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Zb  FORHAK  2(  Fib.  5.5X  ) , 3(  F20.6,5XM 

26  FORMArdH  , //AX.  16HRA0IUS  ( Ml  CRONS  I . AX  , 1 AHS  i Z RAKAMETER 
«Exn  ,17X.7HO  (SCA),  19Xt9HU  IRAnARll 

27  FORMAT (/*  NORMALIZATION  FACTOR  FOR  SIZE  DISTRIBUTION  * • 

28  FORMAT!/*  KIEXT)  = A.EIB.T,*  K(SCA)  =*,£13.7,*  X(RAD) 

• /) 

29  FORMAT!/*  WAVENUMBER  = *,£16.7,*  CM- 1*  ,SX , aDENS I T ¥ = *,£ 
*TICLES  PER  CM-3*/l 

END  PRELIMINARIES 

IFKKK  .GT.  1)  GO  TO  S2 
READ! b, b) WAVE, OFNS 
IF|UENS.EiJ.O.O)l)ENS=l.U 
GNU=1 .OE+OA/WAVE 
Nr ANGE= 1 

REAU(b, 1) lUSTP.NRAOI .NCROS 
I FI IUSTP.EO.b)NRANGE  = NRADI 
ITMM=I T-1 

PI=3.1^1b926b35898 
WRIT  El  6, 15) 

IFI IDSTP.NE.O IGO  TOI 3A,  36, 3 8,AO,A3,A7, AO,51  ) , lOSTP 

READ! b, A) ( F(J ) ,R( J ) , J=1,NRADI 1 

WRITE(6,21I 

WRirE(6,3l(R( J) ,F(  J),  J=l,NRAOI ) 

FSUM=0.0 

00  33  J*1,NRA01 

33  FSUM.FSUM+F(  J) 

GO  TO  49 

34  READ! b,4|RBAR, SIGMA 
WRirE(6,16IRBAR, SIGMA 
RHI=RBAR*  EXP(3.0*SIGMA) 

RLO=RBAR*  EXPI -3.0*S1GMA) 

RADS=  FLOATINRADI-l) 

OELROr(RHI-RLO)/RADS 

FSUM=0.0 

0EN=2.0*SIGMA* SIGMA 
00  35  J=1,NRA0I 
RJ=J-1 

RIJ)=RL0+RJ*nELR0 
GNUM=A  log  I R IJ ) /R  BAR ) 

FIJ)»  EXP(-GNUM*GNUM/OENI*RBAR/R( J) 

35  FSUM=FSUM+FIJ ) 

GO  TO  49 

36  REAUI b, 5IRL0,PHI ,CUE, A,B 
WRITE! 6, 17)RL0,RHI .CUE, A,B 
RAOS=  FLOATINRAOI-n 
0ELR0=!RHI-RL0)/RADS 
FSUMaO.O 

00  37  J=1,NRA0I 
RJ-J-l 

RIJ)*RLO+RJ*DELRD 

F!J)*CUE*A*  EXP! -A*R! J ) I ♦! 1,0  -CUE)*B*  EXP!-B*R!JII 


00  22  SO 
002260 
,13X,7HO  I0022Z0 
002280 
,E1A.7/I  002290 

» »,E13. 7002300 
002310 
15.7,  *PAR002320 
002330 
002340 
002350 
002360 
0023ZO 
0023R0 
002390 
002400 
002A10 
002A2() 
002430 
002AA0 
002AS0 
002460 
002A70 
002AB0 
002A90 
002500 
002510 
002520 
002530 
002540 
002550 
002560 
002570 
002580 
002590 
002600 
002610 
002620 
002630 
0026A0 
002650 
002660 
002670 
00  2 680 
002690 
002700 
002710 
002720 
002730 
0027AO 
002750 
002760 
002770 
002780 


i 


I 


37  FSUM-FSUK+Fl  J) 

002790 

GO  TO  49 

002800 

38  FSUM«0.0 

002810 

0ENS>1.378E>04 

002820 

0ELR0«0.02 

002830 

DO  39  J-l,NKAOI 

002840 

RJ-J-1 

002850 

R(J)«0.02  4-RJ«OE>.KO 

002860 

IF(  J.LT.S)F(J>ii4S0.2 

002870 

IF(J.GE.S)F( J)32.251  *0ELKD*R ( J)««(-4.0 ) 

002880 

39  FSUM-FSUH+F(  Jl 

002890 

GO  TO  49 

002900 

40  IF( IDSTF.EU.4)  GO  TO  41 

002910 

KEAD(!>,5)  VIS 

002920 

RLn>0. 1 

r 

002930 

KHI«1!>.0 

002940 

CUE=30.0 

002950 

A»4.0 

002960 

NRA0I>300 

002970 

DENS>1.1  *10.0  **(+5.0  -ALOGIOIVIS  ) ) 

002980 

41  IF( IUSTP.E0.4)  «EA0(5,4)  KLO.RHI .CUE, A 

002990 

WRtTE(6,18)RL0,KHI,CUE.A,VIS 

003000 

FSUfUU.O 

003010 

«AOs«  float{n«aoi-h 

003020 

0ELR0=<RH(-RL0)/KA0S 

003030 

0042  J*1.NRADI 

003040 

RJ-J-1 

003050 

mj)«RLO+i<J*D£L«0 

003060 

F(  JlaCUE*^  J)**(-A) 

003070 

IF( I0STP.E0.7.AN0. J.LE.3)  F( Jl= l.OE+05 

003080 

42  FSUM=FSUM+F(  J) 

003090 

GO  TO  49 

003100 

43  JZ=0 

003110 

FSUMsU.O 

003120 

00  45  ISZ=1,NRANGE 

003130 

REAO( 3, 5)RL0,RHI .KC.ALF.GAM.OENS.NRAOI 

003140 

WRITE! 6, 19)RL0,RHI ,RC, ALF, GAM. NRAOl , JZ 

003150 

FSUMniSZ)»0.0 

003160 

DELR0>(RHI-RL0)/  FLOAT(NKADI-l) 

003170 

B=ALF/(GAM*«C**GAM 1 

003180 

00  44  J=l,NRAOI 

003190 

R«  J+JZ)=-<LO+OELRO*  FLOAKJ-l) 

003200 

IF( J.EO. 1)XMN=2.0  *PI*K( J+JZ)/(WAVE) 

003210 

I F( J. EU.NRADl )XMX  = 2.0  *PI *R ( J*JZ ) / ( WAVE  ) 

003220 

F( J+JZ )=(  exP(-B*«( J+JZ I**GAMI*R( J+JZ)**Ai,F»*OELRD 

003230 

44  FSUMniSZ»»FSUMI(ISZ)+F(J  + JZ) 

003240 

DELR0bDELK0*2.0  *PI/(WAVEI 

003250 

WRirE(6,24)XMN,0ELR0.XMX 

003260 

45  JZ-NRAUI'fJZ 

003270 

DO  46  JNbI.NRANGE 

003280 

46  FSUM«FSUM+FSUMI  ( JN) 

003290 

NRADI>JZ 

003300 

GO  TO  49 

003310 

47  READ(5,4)  VIS 

003320 

I 


3A 


FSUM.0.0 
0ELLKs0.02 
A*  500. 

NRA0I=280 

GAM»3.0 

HR  I TE( 6. 20) A, NR ADI , GAM, VIS 
PHR=-2.0  •GAM/3.0 

R(l)=  (VIS-0. 04  )/0.31 

F(  1)  = A*R(1)**PW.< 

0048  I=2,NRADI 
EXP0=AL0G(R( I-l) )+0ELLR 
R(I)=  EXP(EXPO) 

F(  I ) = A*«  ( I 

48  FSUM=FSUh+F( 1 ) 

49  On  So  J=1,NKA0I 

50  F( J)=F(J)/FSUM 
WRITE) 6,271 FSUM 
GO  TO  S2 

51  HRITE(6,7) 

GO  TO  61 

52  0TSUM=0.0 
OSSUM=0.0 
0RSUM=0.0 

DO  53  J=1,IT 

53  PSUM(J)=0.0 

IF(<KK  .or.  1)  GO  TO  100 
REAO(  5,  DNINOX 
DO  57  NK»1,NIN0X 
REAO( 5,2)EM,CAY,EMM,CnNC 
CONC  = 1.0 

WRITE! 6,22) NK, EM, CAY, EMM,CONC 
CAY=CAY/eM 

I F(MORTE.EO. 12345) WRITE  (6,26) 

100  CONTINUE 

00  56  L=1,NRA0I 

ALPHA=2.0  MPI^EMM  YRII-I/WAVE 

CALL  FMI EG2( EM, CAY , ALPHA ,0T , OS, OR , P, MCRTE ) 
IF(M0RT£.E0. 1234  5) W« I 76(6,25 )R ( L ) , ALPHA ,0T ,0S , OR 
DO  54  J=1,IT 

PSUM( J)3PSUM( J) +P( J)»F(L)*C0NC 

54  CONT  INUE 

55  OT  5UM=0 TS  UM+R ( L ) ‘R ( L ) *0T4F ( L ) •CONC 
0RSUM=0RSUM+R(L)*R(L)*0R*F(  L)*C0wC 

56  OS  SUM=OSS  UM+R ( L ) *8 ( L ) *0S*  F ( L ) •CONC 

57  CONTINUE 
OTSUM=OTSUM*PI 
ORSUM=ORSUM*PI 
OSSUM=OSSUM*PI 
alb D0=0SSUM/0TSUM 
CAYNG=-CAY 

PFACT=  WAVE*HAVE/(PI*0TSUM#EMM*EMM) 

0058  J=l, IT 

58  PSUM(J)s  PSUM( J)*PFACT/( 4.*PI*AL«00) 

I F(NIN0X.GE.2)HR1T6(6,23) 


003330 

003340 

003350 

003360 

003370 

0033R0 

003390 

003400 

003410 

003420 

003430 

003440 

003450 

003460 

0034T0 

0034R0 

003490 

003500 

003510 

003520 

003530 

003540 

003550 

003560 

003570 

003550 

003600 

003610 

003620 

003630 

003640 

003650 

003660 

003670 

0036B0 

003690 

003700 

003710 

003720 

003730 

003740 

003750 

003760 

003770 

003750 

003790 

003800 

003510 

003820 

003830 

003840 

003850 
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WRIT  E(  6«8  ) I DSTPi  WAVEi  EM  , CA  YMR,  OrSItM.OSSlJM,  ALBOO 
OT SUM=QTSUM*0ENS» 1 .OE-06 
OSSUMsOSSUM«OENS«1.0E-OA 
OR  SUM=ORS  UM*OENS* 1 . Oc -06 
WRITE( 6.2B)0TSUM.0SSUM,0RSUM 
MRirE(6,29)GNU,0ENS 
HRITE(6,9) 

WRITE(6, 10)( (C( J),PSUM( J) IT) 

IF  (lOUT.LT.nr.O  TO  61 
DO  59  J*1,I0UT,4 

WRITE!  6>  10)C(  J)»PSl)M(  J)  ,C(  J + 1)  fPSIIMI  J+1 ) »C  ( J + 2 ) . PSUM  ( J+2  ) ,C  ( J+3)  , 
*PSUM( J+3» 

59  CONTINUE 
return 

61  STOP 
END 

SWROUTINE  FMIEG2(  EM,  CAY  , A'.  PHA  , SOT,  SGS  t SGR  , P ,m:  RTE  t KKK  ) 
dimension  C( 100),EYE1( 100),cYE2( 100),P( 100) 

DIMENSION  REAN(250) ,REBN( 25U) ,FAN( 250) ,FBN(250) 

C0MM0N/BK2  / C,  ALBOD, LULL ,NCROS , I T , I TT ,NR AD  I 

P I E=3.  I'H  59265358  98 
EN-1. 

S = l. 

ISW1>1 
SUMT=0. 

SUMS»0. 

SUMRR.O. 

SUMRI-0. 

SUMSl-0. 

SUMS2«0. 

C MIE  series  cutoff  CRITERION 
FACT=1.2 

IF( alpha. GT. 51. 0)FACT=1.0  +2.26  *ALPHA*A( - .6 13 ) 

A«EM«ALPHA 

B>A«CAV 

gamma-em*cay 

SINA=  SIN(A) 

C0SA>  COS(A) 

COSHB=(  EXP(B)+  EXP(-B) )/2. 

SINHB=<  EXP(B)-  EXP(-B))/2. 

AB«A*A*B*8 
RNLlsSI NA«C0SHB 
SNL1>-C0SA«SINHB 
TNLl=  SIN(ALPHA) 

UNUl=  COS(ALPHA) 

C THESE  ARE  THE  BESSEL  FUNCTONS  OF  N» 1 ORDER 
RNcIRNLlOA-SNLl^BI/AB-COSA^COSHB 
SN>(RNLl«B'»SNi.l*A)/AB-SINA»SINHB 
TN-TNLl/ALPHA-UNLl 
UN-UNLl/ALPHA-i-TNLl 
GO  TO  2 

1 TWONLl-2.0  ♦EN-1. 

STEP7«TW0NL1/AB 

STEP8-TW0NL1/ALPHA 


003860 
003870 
003HHO 
003890 
003900 
003910 
003920 
003930 
003940 
003950 
003960 
003970 
003900 
003990 
0040f ') 
004010 
004020 
004030 
004040 
004050 
004060 
004070 
004080 
004090 
004100 
004110 
004120 
004130 
004140 
004150 
004160 
004170 
004180 
004190 
004200 
004210 
004220 
004230 
004240 
004250 
004260 
0042  70 
004280 
004290 
004300 
004310 
004320 
004330 
004340 
004350 
004360 
004370 
004380 
004390 
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THESE  ARE  THE  RECURRENCE  RELATIONS  EOR  THE  IINRRIMEO  FUNCTIONS  ANO...  004400 

RN=STEP7*( A*RNLI-B*SNL1)-RNL2  004410 

SN=STEP7*(  B*RNLl+A*SNLl»-SNLiJ  i)044/>0 

TN=SrfcP8*TNi.l-TNL2  004430 

UN=STEP8*UNL1-UNL2  004440 

...FOR  THE  PRIMED  FUNCTIONS  004450 

2 RPN=RNL1-EN*( A^RN-B^SNI/AB  004440 

SPN=SNL1-EN*(B*RN+A*SN)/AB  004470 

TPN=-EN*TN/ALPHA-»-TNLl  0044R0 

UPN=-EN*UN/ALPHA+UNL1  004490 

REFERENCE  LIGHT  SCA  TTER  I N^i  . . . BY  VAN  DE  HULSf  004500 

P123  FOR  AN  ANO  BN  004510 

N=EN  004520 

E< SN=KN»TPN-EM*RPN»TN-GAMMA*SPN» TN  004530 

CAPN=SN*TPN-EM*SPN»TN+GAMMA*KPN*fN  ' 004540 

EP  SN=KN4TPN-SN*IIPN+EM#  I SPN»UN-RPN*>TN)  -GAMMAS  (RPN*UN+SPN*TN  ) 004550 

PH  I NsRNOUpN  + SNSTPN-EM# I KPN»UN+SPN»  TN) -GAMMA* I SPN*0N-RPN*TN ) 004  540 

E<  SPN=RPNSTN-bM»RNaTPN-GAMMA»SNSTPN  0045/0 

CAPPN=SPN*TN-EM»SN*TPN+GAMMA»RN*f PN  004580 

EPSPNsRPNSTN-SPNaUN-cMSIRNSTPN-SNSUPN ) -GAMMA*! RN*UPN+SN*T PN ) 004590 

PHIPN=RPN*UN+SPN*TN-EM*(  RN»UPN  + SN»rPN)  •t-GAMMA*(RN*TPN-SN*UPN)  004400 

W=Ei<SN/PHIN  004410 

these  keep  the  ANOS  ANO  BNOS  CLOSE  TO  1, SINCE  THE  NUMERATORS  ANO  004420 

DENOMINATORS  ARE  INOIVIDUALLV  LARGE  004630 

X=CAPN/EPSN  004440 

Y=CAPN/PHIN  004450 

2=EKSN/EPSN  004660 

OENOM=EPSN/PHIN+PHIN/EPSN  004670 

WP=EKSPN/PHIPN  0046RO 

XP=CAPPN/EPSPN  004690 

YP=CAPPN/PHIPN  004700 

ZP=EKSPN/EPSPN  ^ 004710 

OENMP=EPSPN/PHIPN+PHIPN/EPSPN  004720 

REANIN)  = IWP+XP)/OENM.'  004730 

REBN(N)=(W+X)/DENOM  004740 

FAN(N)=IYP-ZPI/OENMP  004750 

FBN(NI=(Y-Z)/DENOM  004740 

I F(MCRTE.E0.67890I WRl TEI6, 19 IREANIN) , FAN! Nl ,REBN(N) ,FBN(N ) 00^.770 

19  format  I 1X,4( E20.  10, lOX ) I 0047R0 

S*-S  004790 

TONP1=2.0  •EN+1.  004000 

SEE  P127  004R10 

SUMT  = SUMT  + T0NP1»(REANI  N) ♦REBNI  Nil  0040 20 

SUMS=SUMS  + rONPl*IREANINI*REAN|NI +F AN I N ) * F AN ( N ) ♦REBN ( N ) *RE BN  I N ) +FB  0040  30 
1N(N)*FBNIN) ) 004040 

SUMRR=SUMRR+T0NPI*S*(REAN(N)-K£6N(N) ) 004050 

SUMRI  = SUMKI+T0NP1*S*(  FAN!  NI-FBNIN) ) > 004040 

IFIEN.EO.l.OIGO  TO  6 004870 

EK=EN-1.  004000 

SUMS1  = SUKS1  + EK*I  EK  + 2.  I * < K AN*  REAM  N) +RBN*KEBN,I  N » +PAN*F4N  !N  > +PBN004890 

**F3N!N))/EN  004900 

SUMS2  = SUMS2  + ! 2.0  *EN-l.  )*  ( RAN*RBN+PAN*PBN'|  / ! EN*EK ) 004910 

6 RAN*REANIN)  004920 

RBN=REBNIN)  004930 
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PAN=FAN(N) 

PBN=FBN(NI 

TERMN=  ABS(REAN(N)I+  ABS(FAN(N|)+  ABS(  REaN(N  ) j'f  ABS(FBN(NI) 
IF(rEKMN.GE.l. 0E-9.AN0.EN.lt. 8.0  ♦FACT*ALPHA )G0  TO  A 
GO  T0(3,t>),ISWl 
3 1SW1>2 
A EN>EN-fl. 

RNL23KNL1 

SNL2“SNL1 

TNL2=TNL1 

UNL2-UNL1 

RNL1=KN 

SNL1=SN 

TNL1  = TN 

UNL1=UN 

GO  TO  1 

5 ALF2=ALPHA*ALPHA 

OEXT*=SGT  OSCAT=SGS  OABS=SGA  ORAnAR  = SGK  AVE  CnS(O) 

0P«=SGMP(SEE  P128) 

^SGT=2.0  •SUMT/ALF2 

SGS>2.0  ♦SUMS/ALF2 

SGA«SGT-SGS 

SGR=  ( SUMRR»SUMR><  -f  SUNK  ( « SUNK  I > /ALF2 

SGMAS=2.0  •( SUMSl+SUMS2)/SUMS 

SGHP=SGT-SGMAS»SGS 

00  98  KsltlT 

COST»C{K» 

SINT*  SORT! 1.0  -COST*COSTI 
PINLI*0. 

RHNLlscO. 

PIN*1. 

RHN>0. 

TAUN=COST 

SUMIlrt=0. 

SUNIlIsO. 

SUMI2R=0. 

SUHI2I*0. 

SUMMATION  OF  MIE  SERIES 
00  97  L=1,N 
EN»  FLOAT! L) 

TM0NL1=2.0  *6N-l. 

ENNP1*EN«( EN+1. I 
ENL1»EN-1. 

I F( L.EU. I )G0  TO  10 

PIN=(T  W0NL1AC0ST’»PINL1-EN*PINL2)/ENLI 
RHN*TW0NL1*PINL1+RHNL2 
TAUN*PIN*COST-SlNT*SINT*RHN 
10  T0NP1=2.0  ‘EN+l. 

SUMI1K  = SUMI  IR  + TONPlAIREANILJ^PIN-t-REBNIDATAUNI/ENNPl 
SUMm*SUMIlI+TONPl*(  FAN(L)*PIN+FBN(L)*TAUN  j/ENNP  1 
SUMI2R«SUMI2R+T0NPI*(REBN(L»*PIN*REAN(L J*TAON l/ENNPl 
SUM12I*SUMI2I+T0f*»l*(  FBN(  DAPIN* FAN(L»*TAUN  »/ENNP  1 
PINL2*PINL1 
PINL1*PIN 
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ooAqAo 
ooAQsn 
00A9A0 
00A97() 
00A930 
00A990 
005000 
005010 
005020 
005030 
00  50  AO 
005050 
005050 
005070 
00  50R0' 
005090 
005100 
Sr,M0051 10 
005120 
005130 
0051A0 
005150 
005150 
005170 
005150 
005190 
005200 
005210 
005220 
005230 
00  52  AO 
005250 
005250 
005270 
005250 
005290 
005300 
005310 
005320 
005330 
0053AO 
005350 
005350 
005370 
005350 
005390 
OO5A0O 
005A10 
005A20 
005A30 
005AAO 
005A50 
005A50 
005A70 


RHNL2=RHNLI 
97  RHNL1=RHN 

11  EYE1(K)=  (SUM!  IR^SIJMI  IR  + SUMt  ir*StlMIH) 
£YE2(  K)=  <SUf1I2R*SUMI2K+SUMI2  l*St)M12I  ) 
P(K»=  ( EYEK Kl ♦EYE2( K1 »/2.0 

9a  CONTINUE 

return 


005AR0 

005490 

005500 

005510 

005520 

005530 

005540 
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